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COMPARISON OF METHODS FOR IDENTIFYING PILOT DESCRIBING FUNCTIONS 


FROM CLOSED- LOOP OPERATING RECORDS 
Rodney C. Wingrove 
Ames Research Center 


SUMMARY 


This report considers the problem of identifying the relationship between 
the pilot's input and output when he performs routing tracking tasks while 
controlling a system with dynamics approximating those of an airplane. A 
difficulty in using only the measured input and output data is that any extra- 
neous output noise by the pilot is transferred through the control loop and 
results in an identification bias error. In this report three different 
identification methods are used to compare the amount of bias error. The 
three methods, which are the parameter model, orthogonal filters, and impulse- 
response techniques, are applied to the identification of both simulated 
(i.e., known) systems and piloted systems. 

The identification of simulated systems illustrates that the amount of 
bias error depends upon the assumed system model. The bias error is reduced 
when the model used in the identification method has the same form (i.e., same 
number of coefficients, dynamic elements, etc.) as that of the simulated 
system. The parameter model method incorporates such a restricted model and 
consistently has the minimum bias error. The orthogonal filters and impulse- 
response methods allow, respectively, more general model forms and have more 
bias error. 

The three identification methods are shown to estimate pilot describing 
functions adequately from representative tracking task data from both single- 
input and two-input control tasks with various levels of external disturbances. 
The results appear satisfactory even when the primary excitation comes from 
the pilot's own output noise. 


INTRODUCTION 


The input-output response of a pilot must be regarded as random, 
nonlinear, and dependent on the task he is performing. Many previous studies 
have shown that this type of response can be represented appropriately wi-th a 
quasilinear system modeled by a linear element (describing function) and a 
remnant term (output noise) . The pilot describing functions usually have been 
identified from records obtained in ground-based simulators (ref. 1) and 
flight tests (ref. 2) wherein carefully controlled external forcing functions 
are used to excite the pilot-vehicle system. The pilot describing functions 
are measured by comparing the input and output signals of the pilot with the 



known forcing function. This technique minimizes those errors in 
identification due to any correlation of the input signal with the pilot’s 
noise.*’ Reference 3 is a good review of this work and summarizes the measured 
pilot describing functions. 

Nfost other methods (refs. 4-12) for measuring pilot describing functions 
depend on random disturbances (e.g,, aerodynamic turbulence, propulsive dis- 
turbance) to excite the pilot-vehicle system. These methods compute the 
describing function of the pilot from only his input and output signals. How- 
ever, there is a fundamental error in identification because the output noise 
of the pilot is transmitted through the control loop and correlates with the 
pilot’s input and output signals, thereby causing a bias error in identifica- 
tion. Reference 4 has shown that this bias error will be small if the ampli- 
tude of the pilot's noise is small compared with the amplitude of the other 
random disturbances in the control loop. 

During routine flight-test operations, there are no carefully controlled 
forcing functions and even the random external disturbance may be quite small 
so that the principal system excitation may come from the pilot’s output noise. 
Reference 5 has shown that, in this situation, it still may be possible, under 
certain conditions, to determine the pilot describing function without incur- 
ring an unacceptable identification error. Reference 5 showed that when the 
computer processing is constrained to identify only physically realizable 
describing functions, the identification error can be reduced by shifting the 
input signal during the computer processing an amount approximately equal to . 
the time delay of the pilot. The theory shows that the identification error 
can be made small, by the time shift procedure, if the correlation time of 
pilot’s noise is small (i.e., near "white" noise) compared with the sum of all 
time delays through the control loop. A somewhat different approach in refer- 
ence 6 indicated that if the form of the identification model (e.g., number of 
coefficients, dynamic elements) used in the computer processing were 
restricted to be the same as that of the actual pilot’s describing function, 
then the identification bias error may also be reduced. The studies in refer- 
ences 5 and 6 have shown, in effect, that the identification error may be 
strongly affected by the constraints in the identification model. 

This report compares the results obtained from three identification 
methods which differ in the constraints placed on the describing function 
model. The impulse response method (refs. 4, 5, 10, and 11) assumes only that 
the describing function is physically realizable (i.e., the impulse response 
is constrained to be zero for negative time, ref. 5). A more constrained 
model is obtained with the orthogonal filter method (refs. 7 and 9) which 
assumes, in addition, that the describing function can be represented by a 
finite series of orthogonal filters. The most constrained model is obtained 
with the parameter model method (refs. 6, 8, 12) which assumes that the 
describing function can be represented by one specified model form (a simple 
second-order model is used for the results in this report) . These three 
methods will be applied with the time shift procedure used in reference 5. 


2 



This report first describes the three identification methods and applies 
them in the identification of both simulated and piloted systems. The identi- 
fication of the known systems will be used to show the effect of the output 
noise (relative magnitude and spectrum) on the identification bias error. 

These results will be compared with the theory of reference 5 and will- also 
be used to show how the constraints, incorporated within each method, a.ffect 
the identification error. The identification of piloted systems will illus- 
trate the application of these methods for representative single-input and 
multi -input tracking tasks. 


IDENTIFICATION METHODS 


Before outlining the three identification methods used in this report, 
we shall briefly discuss the piloted control system elements. Figure 1 
presents a block diagram of the pilot in a representative compensatory 
tracking task trying to control his output y(t) so that the input error 
signal x(t) is kept near zero. Generally, the input-output characteristics 
of the pilot must be considered as complex, nonlinear, and time varying. How- 
ever, for the purposes of modeling, it is common practice to assume that his 
characteristics can be represented by a quasilinear system (ref. 3), This 
mathematical model contains the linear element G and the noise source n. 

The element G(jo)), which is called the pilot describing function,^ is the 
frequency response of a linear constant coefficient system to the input x(t). 
The term n(t) represents the difference between output of the pilot, y(t), 
and the output of the describing function G(jo)) driven by x(t). Thus, 
n(t) accounts for remnant terms such as nonlinearities, time variations, and 
additive noise in the output of the pilot. 

The controlled system is mathematically characterized by the constant 
linear element H and the noise source i. The time history i(t) accounts 
for nonlinearities and time variations in the controlled element, time-varying 
commands, and all disturbances from aerodynamics, propulsion, etc., external 
to the pilot. 

In this report, we shall use each of^the three methods to compute from 
the records x(t) and y(t), an estimate GCjus) that represents the best 
linear relationship between x(t) and y(t). "Best" here means that the 
integral of the squared residual. 


^Technically, G(jco) represents a random input describing function 
because random, rather than sinusoidal, signals are used here (see ref, 3) , 
Also, to avoid additional notation, terms such as G(joj) and H(jto) will be 
used to represent both the transfer functions of linear systems and the 
describing functions of nonlinear systems. 
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is minimized over a given record length, where e(t) is the difference between 
the actual record y(t) and the output of the system G(jto) excited by x(t) . 
The fpl lowing discussion will briefly describe the formulation of the three 
methods for the identification of single-input single-output systems. The 
formulation for multi-input systems (e.g., two input, one output) along with 
the details of the computer programing is further described in appendix A. 


Parameter Model 

The parameter model method assumes a particular describing function model 
for the pilot dynamics and then solves for the parameter in that model. The 
model^ used for the results in this report has the form 


GCjto) 


a3ju) + -Ajo) 

(ju ))2 + aijo) + &2 


( 1 ) 


Estimates of the parameters ai, Si2> ^ 3 > determined, as shown in 

appendix A, by a quasilinearization technique. The time shift X accounts 
for any pure time delay in G(jo)). A method for choosing an appropriate value 
for X is illustrated in appendix B. 

The parameter model method is restricted in that only a limited set of 
systems, which have the specified form of equation ( 1 ), can be adequately 
identified. 


Orthogonal Filters 


The orthogonal filter method is somewhat more general than the parameter 
model method. It assumes that the unknown system dynamics can be modeled by 
a series of transfer functions of the form 


G(jm) 


= 


bi 

Tj jco+l 


baCTijm-l) 
(Tijm+l) (T2jto+l) * 


bsCtijoj-l) (T 2 jo)-l) 

(Tljto+l) (T2jo)+l) CTajw+l) * ' J 


( 2 ) 


This series is commonly called a set of Kautz (or linearly independent) 
filters. Estimates of the parameters bi, b 2 , bs,. . . etc., are determined, 
as shown in appendix A, by a multiregression technique. For the results in 
this report, the first five filters in the series were used. The values of 
the time constants T 2 , T 2 , . . . T 5 were taken (following ref. 7 ) as 1 , 0 . 5 , 
0.25, 0.125, 0.0625 sec, respectively. 


^This model form appears to be reasonable for the results in this 
report; however, the best model form will depend somewhat on each particular 
situation. 
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Impulse Response 


The impulse response method assimes a very general input-output 
relationship that can be represented by the form 

GCju) = f g(T)e~^^“ dx (3) 


Here g(x) is an impulse response function that is assumed to be zero for 
X < 0 (i.e., g(x) is constrained to be a physically realizable system) and 

also zero for t > Xj„ (i.e., a finite memory time Xj„) . The estimates for 
the impulse response function, as described in appendix A, are calculated at 
discrete times, g(0), g(At), g(2At), etc. For the results in this report, the 
estimates were made at 20 discrete points. 


IDENTIFICATION OF A KNOWN SYSTEM 


In order to illustrate the results to be expected from each method just 
discussed, we shall first consider an example where the pilot dynamics are 
simulated; that is, where the describing function to be identified is known. 
For this example, the elements G and H were^ 


and 


GCjco) 


40 + 50 

(ju)2 + 3.5 jo3 + 25 ® 


HCjo)) 


j cj + 0 . 585 

jo) TO (jj) ^ + 0 . 592 jw + 0 . 584 J 


The dynamics for this example were simulated on a digital computer. The 
output of a random noise program was appropriately filtered to obtain the 
desired spectrums for n(t) and i(t). The resulting records for x(t) and 
y(t) were used with the^three methods, as described in appendix A, to 
determine the estimate, G(jm). 

Figure 2 presents the computed magnitude |G(jco) | and phase angle 
<); GCjoj) as functions of frequency. Also shown for comparison are the magni- 
tude |G(ju) I and phase angle G(jm) of the actual system. Two cases were 
simulated in order to illustrate the effect of the excitation source on*the 
identification error. In the first case, the system dynamics are excited by 
both the external noise source i(t) and the internal noise source n(t). 

In the second case, the dynamics are only excited by the internal noise source 
n(t) and i(t) = 0. 

^The dynamics of G and H were taken to be near those of the piloted 
system in the next exairple. 
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Figure 2(a) presents the results for the case in which the system is 
excited by both the internal noise source and the external noise source. This 
figuj-e shows that with this moderate amount of external disturbance 

(i^ = 0.5 n^) the estimated describing functions GCjm) are generally near 
the actual system describing function The estimate derived by the 

paraiheter model method provides a nearly perfect match with the actual system 
G(jo)). The estimates derived by the orthogonal filters and impulse-response 
methods show only small differences between SCjco) and GCju). These differ- 
ences are probably due to modeling error these methods cannot exactly 

match the actual system model), and, because there is a bias error, due to the 
feedback of n(t)^ through the control loop. The effect of this bias error 
is discussed in more detail with the next case. 

Figure 2Cb) presents the results for the case in which the system is 
excited by^ the internal noise source n(t) and there is no external distur- 
bance, i(t) = 0. This figure shows that with only an internal excitation, 
GCjm) may be quite different from GCjto). This difference, or ''identifica- 
tion bias error," is due to (ref. 5) the noise n(t) being transmitted 
through the control loop producing a correlation between n(t) and x(t), so 
that the estimate tends to measure the negative inverse of the alternate path, 
-1/H, rather than measure the actual system G. The identification bias 
error is most clearly seen for the measured magnitude shown on the top of 
figure 2(b). The estimates |G(jo)) | generally tend away from the actual 
system jG(ja))| toward the curve shown for |l/H(jo))j. 

The analysis in reference 5 has shown that the amount of bias error 
depends upon the relationship between the correlation time, of the noise 
and the time delay, xg, in the control loop. An analytical expression for 
this bias error can be written as 


G(jo)) 


G(jm) 


+ e 

V 


-W^G 


G(jto) - ^(jw)_ 


•TG(j‘^) 


actual system 


bias error 


( 4 ) 


The primary assumptions in formulating this expression (ref. 5) are; 

(1) i(t) = 0, (2) X = Tg, (3) the estimate G(jo)) is constrained to be 
physically realizable, and (4) H(ju)) is minimum phase. 

The magnitude of the bias error is directly related to the constant 

factor por the simulated examples shown in figure 2, both the cor- 

relation time Tjj of the noise and the time delay xg in the simulated 
system were taken as 0.2 second. In order to illustrate the effect of these 
two ^'quantities on the identification error, several runs were made with other 
values of x^^ and Xg. Figure 3(a) presents results for various values of Xjj 
with Xg held constant at 0.2 sec. Figure 3(b) presents results for various 
values of xg with Xj^ held constant at 0.2 sec. The estimated magnitude 
1 g| is presented at one frequency, w = 1 rad/sec. These estimates from the 
three identification methods are shown in comparison with the theory (eq.(4)) 
and with the actual value for the system (|g| = 8.4 dB) at w = 1 rad/sec. 
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The relative effects of noise correlation time Tj^ and the time delay 
tq are apparent in figure 3, Increases in the noise time constant are 

seen in figure 3(a) , to increase the bias error. The error will only be* neg- 
ligible when Tj^ is near zero, that is, when the noise spectrum is near 
"white." Increasing the time delay tq (fig. 3(b3) decreases the ideptifica- 
tion error. The trends in figures 3CaJ and 3(b) follow the theoretical con- 
siderations that the identification error will be small when the noise time 
constant is small compared with the time delay through the control loop. 

Figure 3 shows that the difference between 1 g| and [g| depends strongly 
on the particular identification method. We can note that the parameter model 
method, in which the model has the same form as the actual system, has an esti- 
mate |g| which is consistently closest to |g|. The orthogonal filters and 
impulse-response methods which are, respectively, less constrained, have esti- 
mates |fi| further from |g|. The curve for the theory (ref. 5) was derived 
for an estimate that is constrained only to being physically realizable. The 
most general identification method, impulse response, is seen to be near this 
theoretical boundary. (The impulse-response method does not match this curve 
exactly because of restricted memory time noted with eq. (3). 3 

The results in figure 3 show that the identification error is related to 
the inherent constraints in the identification model. The identification 
error will be reduced when the model allowed by the identification method is 
restricted to a form that is near that of the pilot dynamics. In many practi- 
cal situations, however, the exact form of the pilot dynamics will be unknown. 
In such cases, the more general methods, such as orthogonal filters and 
impulse response can be used to gain insight and possibly help to decide on a 
reasonable form for the pilot dynamics. The following exan 5 )les illustrate 
applications wherein the combined use of all three methods aid in 
identification . 


IDENTIFICATION OF PILOTED TRACKING DATA 


Recorded data from the simulation study described in reference 13 will be 
used to illustrate the identification of pilot describing functions. In this 
example, the pilot controls the longitudinal axis of a simulated jet transport 
aircraft. Two cases will be analyzed. The first case is a single-input 
tracking task wherein the pilot controls the pitch attitude of the aircraft. 
The second case is a two- input tracking task wherein the pilot controls both 
the pitch attitude and altitude deviations of the aircraft. 

The results to be presented in this report represent the average values 
computed from 12 minutes of tracking data (for subject A in ref. 133. Where- 
ever possible, a comparison will be made between the results computed for this 
paper, and the results computed from the same data in reference 13. 
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Single- Input Task 


I.n a single-input compensatory tracking task the pilot was trying to 
control his manipulator output y(t) so that the attitude error signal x(t) 
(displayed on an oscilloscope) was kept near zero. The simulated aircraft 
dynamic^ H had the form 


HCjo)) 


jco + 0.585 

jo) [Cjw)^ + 0.592 jm + 0.584] 


and i(t) was a sxqjerposition of eight continuous sine waves (table 1) which 
gave a random- appearing signal. 

Pilot describing functions computed from the experimental records are 
presented in figure 4. Curves are shown for the three computational methods 
in this report and for comparison, the results determined in reference 13 are 
also shown. 

Reference 13 used the standard cross-spectral method. The cross-spectra 

(jw) and were computed at each frequency contained in i(t) and 

the pilot describing function was determined by the ratio $. (jo))/$. (joj) . 

xy X X 

This method, in effect, coirelates the input x(t) and output y(t) with the 
known forcing function i(t) and thereby eliminates the bias error in identi- 
fication. The cross -spectral measurements, shown in figure 4, were made for 
three separate run lengths of four minutes each. The cross-spectral ratios 
are seen to have significant scatter in phase angle at low frequencies. (See 
ref. 14 for a discussion of this measurement problem.) Except for this uncer- 
tainty at the lower frequencies, there is generally good agreement between the 
cross-spectral ratio (jw)/4>j^^(jo)) and the estimates G(juj) determined by 

the three methods in this report. 

The residual terms calculated by the three identification methods for 
this example have the following values: 

£^/y^ 

Parameter model 0.145 

Orthogonal filters .146 

Impulse response .144 

The residual terms, shown normalized with respect to y^, represent about 
14-1/2 percent of the pilots output. It is interesting that the residual for 
the parameter model method is about the same as that for the other less 
restricted methods. As further discussed in appendix B, this good agreement 
indicates that there is negligible modeling error with the parameter model 
estimate. Because there is good agreement between the curves of G(jo)) in 
figure 4, and also good agreement between the residual values, it appears that 
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the chosen form of parameter model satisfactorily represents the pilot 
describing function for this example. 

One way of verifying the results for this example is to consider the 
response of the closed-loop system. Let us take, for convenience, the param- 
eter model estimate from figure 4 for the pilot describing function 


G(jw) 


49 ,2 jco + 62 .4 -0 . 2 joj 

Cja))2 + 3.65 jo) + 28.4 ® 


With this function we can then calculate the transfer function between the 
input signal i(t) and the output state, r(t), as G(jaj)HCjto)/ [l+G Cjo))H(jaj)J . 
This curve is coir 5 )ared in figure 5 with the measured cross-spectral ratio 
$ir(jw)/$ii (juj) . Within the scatter of these cross -spectral measurements, we 

can see good agreement between the calculated closed-loop response, using 
G(jo)), and the ratio (joi) . 


Two- Input Task 

Figure 6 illustrates the two-input single-output tracking task analyzed 
in this case. The pilot is to keep both the inputs Xj (t) and X2(t) near 
zero by the use of the manipulator output yC^). The simulated inner-loop 
element Hj and the disturbance ii(t) are identical to those in the single- 
input task just considered. The simulated outer-loop element H 2 had the form 

juCjw + 0.585) 


and the input disturbance iaCt) was a superposition of seven sine waves 
(table 2 ) at frequencies spaced between those of the eight sine wave 
frequencies in ii(t). 

Estimates for the pilot describing functions that were confuted from the 
experimental records are presented in figure 7. The estimates Gi(ju) of the 
pilot response to the inner- loop signal are presented in figure 7(a) and the 
estimates G 2 (jtii) of the pilot response to the outer-loop signal are pre- 
sented in figure 7(b). These results represent the average values derived 
from twelve minutes of pilot tracking data (again, subject A in ref. 13) . 

In figure 7(a), the cross-spectral results from reference 13 are compared 
with the results computed for this report. Significant scatter can be seen in 
the cross-spectra results. Reference 13 points out that fundamental limita- 
tions restrict the accuracy of these cross-spectral measurements in multiloop 
tasks. The three methods discussed in this report are seen to agree generally 
with each other for frequencies up to about 4 rad/sec. The only major differ- 
ences are in the estimated magnitude |Gi(jto)| at the higher frequencies. 
Generally, in the mid- frequency region between 0.4 and 4 rad/sec, the magnitude 
determined by the three methods agrees with the cross-spectral measurements. 
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and, except for the scatter in the cross -spectral measurements of phase angle 
at the lower frequencies, agrees with the phase angle results. 

The estimates made by the three methods in this report for the pilot 
describing function G 2 (jw) in the outer loop, figure 7(b), show good agree- 
ment. A measurement of G 2 , which might allow comparison with these results, 
was not computed in reference 13. We can, however, check the validity of the 
results in figure 7(b) by considering the estimated dynamics of the closed- 
loop system. Let us take the parameter model estimates (from fig. 7) for the 
pilot describing functions 


f, 14.9 ICO + 14.4 -0.2j 

1(1'^) (j(o)2 + 2.5 jo) + 12.1 ® 

f. ^ 6.3 io) + 3.8 -0.2j 

" (jw)2 + 2.5 jo) + 12.1 ® 


With these values, we can then calculate the estimated transfer function 
between the input signal i 2 (t) and the output state r 2 (t) as 

G2(j^)Hi(joj)H2(jgj) 

1 + Gi (joj)Hi (joj) + G2(jto)Hi(jo))H2(jw) 


This curve is compared in figure 8(a) with the cross-spectral ratio 

(ji»j)/$. . (jto). Although there was some scatter in these cross-spectral 
1 2 ^2 ^2^2 

measurements, we can see good agreement between the calculated closed-loop 
response, when Gi(jto) and G 2 (jco) are used and the computed cross-spectral 
ratio. 

Two external disturbances, ii(t) and i 2 (t) were used to excite the 
dynamics for the two-input tracking data discussed to this point. We will 
next illustrate some identification results when the excitation ii(t) is 
removed. For this situation, which was also simulated in reference 13, esti- 
mates were made with the three identification methods of this report. The 
results for the orthogonal filters and impulse- response methods were found to 
agree well with the parameter model results. Again, for convenience, the 
following estimates derived by the parameter model method only are given 


Gi(jw) 
G2 (jw) 


12.5 jo) + 7.7 -0.2jo) 

(jw)2 + 3.1 ju) + 10 ® 

5.05 jw + 1.8 -0.2jw 

(jo))2 + 3.1 jo) + 10 ® 


The total closed- loop response, calculated from these estimates, is compared 
in figure 8(b) with the corresponding cross-spectral ratio i 
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This figure shows very good agreement between the calculated closed- loop 
response, when GiCjw) and G 2 (ja)) are used, and the 'measured cross -spectral 
ratio. This figure indicates that the results are adequate when only an outer 
loop disturbance, i 2 (t), is used to excite the closed-loop dynamics. 

The identification for a pilot control task is illustrated next for the 
case when there are no significant external disturbances, ij or i 2 - This 
illustration is interesting because only his output noise is used for system 
excitation. 


APPLICATION fo A PILOT TASK WITH NO EXTERNAL DISTURBANCE 


This example uses recorded data from the simulation study described in 
reference 15. This simulation represents a formation flying task in which the 
pilot is trying to hold his aircraft in a constant position behind a lead air- 
craft. Visual cues were presented to the pilot through a closed-circuit tele- 
vision system, and motion cues were applied with the Ames six-degrees-of- 
freedom motion device. For this example, only the lateral -directional 
control task will be analyzed for which the aircraft frequency response can be 
approximated as 


Hi (jw) 
H2(jco) 


4.3(jgj)^ + 2.2 jgj + 6.5 

+ 1.3Cjto)2 + 2 jo) + 2.1J 

0.56 

Cjo3)2 


The function HiCjw) approximates the dynamics from the pilot's lateral stick 
deflection to the aircraft roll attitude. The function H 2 (joi)) approximates 
the dynamics between roll attitude and the aircraft lateral position. In this 
simulation, there were no significant‘s external noise sources, iiCt) or i 2 (t), 
so that the primary excitation came from the pilot's output noise. 

The estimation of two independent functions GiCjw) and G 2 Cjw) presents a 
fundamental problem in this situation because when i 2 (t) = 0, the two input 
signals Xi (t) and X 2 (t) C^'oll angle and lateral position) are exactly corre- 
lated. Solutions with the general identification methods, orthogonal filters 
and impulse response, are not possible in this situation (i.e., the matrix to 
be inverted, eq. (A2) , becomes singular). With the parameter model method, 
however, a solution is possible because the restricted parameters in the model 


“^There may have been a small disturbance from several sources (simulation 
noise, coupling from other axes, long-term drift, etc.); however, the amplitude 
of the disturbance was quite small compared with the pilot's output noise. 
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are still linearly independent (i.e., the matrix, eq. (A5) , is not singular). 
The measurements made with the parameter model give the following results: 


V 3.2 id) + 5.7 - 0 . 2 j 

Cjw )2 + 2.2 jo) + 13.8 ® 

4.3 jo) + 1 .4 - 0 . 2 j 

Cjo))2 + 2.2 jto + 13.8 ® 


GzCju)) 


These results qannot be verified as done in the previous example because 
there are no external disturbance sources, iiCt) or i 2 Ct)» that would allow 
alternate cross -spectral measiixements . It is possible, however, to gain some 
insight into the adeqiiacy of these results from the calculated closed-loop 
transfer functions. These fimctions, calculated with the estimates GiCjw) 
and 62 ( 3 f'j)* presented in the upper portion of figure 9. The transfer 

functions for n to y, n to xi and n to X 2 are presented in figures 9(a), 
9(b), and 9^c), respectively. The calculated closed-loop response, when 
di(jo)) and 62 (jto) are used, shows amplitude peaks at frequencies of 0.5 and 
1.75 rad/sec. These two well-defined peaks correspond to the long and short- 
period closed-loop response of the combined pilot/vehicle dynamics. One would 
expect, therefore, to see large amplitude oscillations at these two 
frequencies in the operating records. 

Power spectra were calculated from the operating records y(t), xiCt) and 
X 2 (t). These power spectra are shown in the lower portions of figures 9(a), 
(b), (c) , respectively. We see that the shape of the power spectra generally 
follow the trends calculated by the closed-loop transfer functions. In par- 
ticular, each shows peaks near the frequencies of 0.5 and 1.75 rad/sec. This 
comparison of the power spectrum with the predicted closed-loop response does 
show that these estimates for GiCjto) and 62 (j^^) t® appear reasonable. 


CONCLUDING REMARKS 


This report has con 5 >ared the use of three methods for the identification 
of pilot describing functions from closed-loop operating records. Each 
method has inherent constraints that restrict the form of the identified 
describing functions. The impulse response method identifies a general class 
of describirig functions that are constrained to be physically realizable. 

The orthogonal filters method identifies a more restricted class which is 
constrained, in addition to the above, to represent a finite series of dynamic 
models. The parameter model method identifies a very restricted class of 
describing functions which are further constrained to represent only one spec- 
ified form of dynamic model. In addition, all of these methods use a 
time shift that is approximately equal to the time delay of the pilot. 

This report shows that the identification bias error due to the 
correlation of the pilot's output noise in the control loop is related to the 
constraints imposed by each identification method. The parameter model method 
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is constrained to a specific form of pilot model and consistently has the 
minimum bias error. The orthogonal filters and impulse response methods are 
respectively less restricted and have respectively more bias error. This 
report also shows that, indepependent of the method, the bias error will be 
negligible only when the correlation time of the pilot's output noise is small 
(i.e., near "white") in relation to the time delay of the pilot. 

The three methods were shown to estimate adequately the pilot describing 
function for representative single-input and two-input tracking task data. 

The primary advantage of these methods (as compared with the standard cross- 
spectral method) is that they can identify the pilot dynamics from routine 
closed-loop operating -t^asks where the external disturbance may be small, that 
is, in situations where the principal excitation comes from the pilot's own 
output noise. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 9403S, Sept. 23, 1970 
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APPENDIX A 


COMPUTER PROCESSING EQUATIONS 


The computer processing equations used for the three identification 
methods in this report are generalized to include any system with several 
inputs xi(t), X£(t) , . . x^^Ct) . . . and with one output y(t). 

Assume that the input and output time histories have been digitized at a 
uniform series of discrete points kAt, where At represents the digitized 
increment. For each of the following techniques the input data are shifted by 
an amount where X represents the time delay of the pilot. The shifted 
input data will be represented as x^'(k), where x^'(k) = x^(kAt - A). This 

time shifting of the digitized data and the processing equations, described 
below, was programmed on a digital computer. The solutions for the impulse 
response and orthogonal filter methods used a multiple regression technique. 
The parameter model method used a quasilinearization technique. 


MULTIPLE REGRESSION 


This technique (ref. 4) minimizes the least squares function 

K 

[yCk) - y(k)]^ 
k=i 



where K is the total number of data points and y(k) is the estimated 
output represented by 


y(k) = 


b.Z.(k) 


j = l 


(Al) 


The terms Zj (k) are the result of operating (as will be shown below) 
on the input data. With this formulation, the coefficient bj can be 
estimated by the following matrix inversion: 
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-1 




K K 


K 

bl 


^ Zi2(k) • • “ ^ ZlCk)Zj(k) 


Zi(k)y(k) 

• 


k=i k^ 


k=i 

• 


K K 


K 



^ Zj(k)ZiCk) • ■ ■ V Z/(K) 


V Zj(k)y(k) 



k=j ■fef 


k-i 


CA2) 


This general computing technique, which involves the inversion of a 
Jxj matrix, was used for the impulse response and orthogonal filters methods. 
These methods differed primarily in the form of the operators, Z., in each 
case. ^ 


Impulse Response 

The operators for the impulse response method are a series of simple 
time delays that operate on the imputs: 

ZiCk) = xl(k) 

ZaCk) = xl(k-i) 

Z 3 (k) = x{(k- 2 ) 


or, in general, 


ZjCk) = xJCk-m) 

where j = (Z-l)M + m + 1, These operators represent a series of M simple 
time delays (m = 0, 1, 2 . . . M-1) operating on the inputs x { (k) , x^Ck). . . 
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x^Ck) .... With these operators, the coefficients resulting from the 
computer solutions (eq. (A2)) are discrete values of the impulse response 
function 


g^Cm) = bj 

where again j = (Z-l)M + m + 1. This time domain solution was transferred 
into the frequency domain by the following sin^ile discrete form approximation 
for the Fourier transform. 


G2_Cju) 


= e 


-Xjcj 


M-1 

m=o 


(m)e 


-mAtjoj 


(A3) 


Orthogonal Filters 

The operators for the orthogonal filters method represent the output of 
a series of filters driven by the inputs: 

M-l 

Zl(k) = di(m)xiCk-m) 

m=o 
M-l 

Z2(k) = d2(m)xiCk-m) 

m=o 

M-l 

Z3(k) = E d3Cm)xl(k-m) 
m=o 


or, in general, 

M-l 

Zj(k) = E di(m)x^(k-m) 
m=o 


where j = (Z-l)I + i and the terms di(ra), d 2 (m) . . , di(m) . . . are a set 
of I discrete impulse response functions representing the set of orthogonal 
filters shown with equation (2) in the main text of this report. 

The coefficients resulting from the computer solution (eq. (A2)) are 
used to calculate the estimate for the system impulse response function 

i^(ro) = Z bjdi(m) 

i=l 
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where j = (Z-l)I + These time domain impulse response functions are 
transformed into the frequency domain by means of equation (A3) . 


QUASILINEARIZATION 


The quasilinearization technique assumes that the pilot dynamics can be 
modeled by a set of constant- coefficient linear differential equations (plus 
the time shift X already mentioned) . If a second-order model is chosen, the 
equations can be written in the form 


— 1 
f — s 

1 

■-ai 1" 

"y(k)" 


3-3 35 . . . 

_w(k)_ 

-a2 0 

_w(k)_ 


.34 3 g . . ._ 


xl(k)“ 
X2 (k) 


(A4) 


where w(k) is a dummy variable and y(k) represents the estimated output of 
the pilot. In order to match the estimated output with the measured output 
K 

(minimize ^ [y(k)-y(k)]^) , an initial estimate is made for ai, a 2 , as- . . . 
k=l 

The following iterative procedure (ref. 16) is then used successively to 
improve the estimate 






K ! K 


K 

3l 


3l 


2?(k) ! zi(k)z2(k) . . . 


zi(k)[y(k)-y(k)] 





k=i 1 k=l 


II 





1 

K ; K 


K 

32 

= 

y\ 

E2 

+ 

D Z2(k)Zl(k); E 2|(k) . . . 


E Z 2 (k)[y(k)-y(k)] 

• 


e 


k=i |k=l 

« « 9 4 9 


» 

« 

new 

9 

^9 ^ 

old 

9 9 « 9 • 


; J 


(AS) 

On each iteration the system differential equations (eq. (A4)) were solved to 
determine the estimated output y(k) and the method of reference 17 was used 
to determine zj (k) , where Zj (k) is the gradient function 9y(k)/8aj. 


For the results in this report, it was found that reasonable convergence 
(from satisfactory initial estimates) was obtained in about five iterations. 
The estimated parameters from this time domain solution can be directly 
related to a frequency domain representation through transfer functions of the 
form 
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Gi (jto) 


asjco + a4 

(joj)^ + aijo) + a2 


(ju)^ + a-ijo) + a2 
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APPENDIX B 


SELECTION OF THE TIME SHIFT A 


The time shift A used with the computer processing equations to account 
for any time delay tq may be approximately known in some situations (e.g., 
ref. 3); but, in general, it will be unknown and will depend on the particular 
piloting task. The purpose of this appendix is to illustrate some results 
from using several values of A and to illustrate the rationale in selecting 
the values of A used for the results presented in the body of the text. 

Figures 10(a) and (b) show the effect of the time shift A on the 

residual (normalized with respect to the output ~y ) . Results for the 

single- input task are presented in figure 10(a) and for the two- input task, 
in figure 10(b). These figures illustrate that the parameter model method is 
more sensitive to A. The orthogonal filters and impulse-response methods, 

which are less restricted, are less sensitive to A. 

With both the single-input and two-input tasks, the residual for the 
parameter model method shows a minimum value at A = 0.2 sec. For this 
method, the value of the pilots time delay tq is one of the parameters that 
must be estimated. It is therefore reasonable to choose a value for 
TG = 0.2 sec. 

The residual error for the orthogonal filter and impulse-response 
methods shows (fig. 10 (a)) a minimum value at A = 0. (This is to be expected 
because with A = 0 there is significant correlation between n(t) and x(t); 
see ref. 5.) The residual is about the same for values of A from 0.1 sec 
to 0.3 sec. It was found that the form of the estimate G(jco), such as shown 
in figure 4, also has no appreciable difference for this range of A from 
0.1 to 0.3 sec. Therefore, for these more general methods, the estimate 
G(joj) is relatively insensitive to the exact value used for A. 

A value of A = 0.2 sec, which is tq estimated by the parameter model 
method, was used for the results presented in the body of the text. 

Figure 10 shows that, near the time shift A = 0.2 sec, all the methods have 
about the same residual (e^/y^ = 0.14 to 0.15) representing about 14 to 15 
percent of the output y^. It is interesting that the amount of residual for 
the parameter model method is about the same as the amount of residual for 
other, less restricted, methods. It appears, therefore, that there is a 
negligible modeling error and the form chosen for the parameter model is 
reasonable for the results in this report. 
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TABLE 1,- THE EIGHT SINE WAVES SUPERIMPOSED TO REPRESENT 
THE SIGNALS i(t) or ii(t) 


Frequency 

Crad/sec) 

0.157 

.288 

.S*!- 

.524 

.969 

1.75 

3.25 

6.00 

11.1 


Relative 

magnitude 

1 

1 

1 

1 

0.1 

.1 

.1 

.1 


TABLE 2.- THE SEVEN SINE WAVES SUPERIMPOSED TO REPRESENT 
THE SIGNAL i2(t) 


Frequency 
(rad/ sec) 

0.209 

.367 

.681 

1.28 

2.38 

4.42 

8.17 


Relative 

magnitude 

1 

1 

1 

0.1 

.1 

.1 

.1 
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(a) Effect of t ; = 0,2 sec. 

n G 


Figure 3,- Effect of x^ and tq on the identification of a known system; 
i = 0, RnnC"'^) "® jA=tq,o)= 1 rad/sec. 
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Frequency, <0, rad/sec 


(a) Two external disturbances; and i^. 

Figure 8.- Estimation of closed-loop dynamics; two-input task 








Phase angle, deg Amplitude, dB 




(b) One external disturbance i^; = 0. 

Figure 8.- Concluded, 
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Frequency, w, rad/sec 


(a) The describing function n to y compared with the power spectrum of y(t) 

Figure 9.- Calculated describing functions compared with measured power 
spectrums; two-input task with no external disturbance. 






Frequency, tj, rad/sec 


(b) The describing function n to compared with the power spectrum of xi(t) 

Figure 9.- Continued. 







Power spectrum Amplitude , dB 




(c) The describing function n to X2 compared with the power spectrum of X2(t). 

Figure 9 .- Concluded. 


3 











National Aeronautics and Space Administration 
Washington, D.C 20546 


official business 


FIRST CLASS MAIL 


postage and fees paid 
national aeronautics ANl 
space administration 


POSTMASTER: 


If Undeliverable (Section 158 
Postal Manual ) Do Not Return 


''The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the aWiosphere and space. The Administration 
shall provide for the xvidest practicable and appropriate dissemination 
of information concerning its activities and the results thereof. 

— National Aeronautics and Space Act of 1958 


NASA SaENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS; ScienriSc and 

technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs^ 
Technology Utilization Reports and 
Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washinston^ D.C. 20546 




